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I ntroduct ion 

The  barotropic  finite  element  code  described  by  Staniforth 
and  Mitchell  in  Ref .  1  has  been  the  point  of  departure  for  exten- 
sive efforts  at  the  Naval  Postgraduate  School  to  develop  improved 
techniques  for  numerical  weather  prediction.  Two  earlier  reports 
(Refs.  2  &  3)  described  applications  of  tensor  product  techniques 
to  extensions  of  the  code  developed  by  Hinsman  (Ref.  4).  The 
present  report  is  concerned  with  proposed  alterations  of  the 
Stani f or th-Mi tche 1 1  code  for  the  purpose  of  improving  computa- 
tional efficiency. 

Three  separate  contributions  are  reported  here.  First  of 
these  is  a  scheme  for  the  direct  solution  of  the  Helmholtz  equa- 
tion as  a  substitute  for  the  Concus  and  Golub  iterative  algorithm 
(Ref.  5)  used  in  the  Stani forth-Mi tche 1 1  code.  The  second  item 
is  an  improved  solution  of  the  generalized  eigenvalue  problem 
that  arises  in  transforming  the  two-dimensional  Poisson  and  Helm- 
holtz problems  into  a  series  of  one-dimensional  problems.  This 
solution  takes  advantage  of  the  special  form  of  the  matrices 
involved  (symmetric,  tridiagonal)  and  avoids  matrix  inversion. 
The  third  item  deals  with  the  special  form  of  the  generalized 
eigenvalue  problem  that  arises  when  there  is  a  periodic  boundary 
condition  in  the  east-west  direction.  In  this  instance  the 
symmetric,  formerly  tridiagonal,  matrices  remain  symmetric,  but 
have  nonzero  elements  in  the  upper  right-hand  and  lower  left-hand 
corners.  A  separate  solution  algorithm  is  required  for  determin- 
ation of  the  eigenvalues  and  also  for  the  eigenvectors. 


Pi  rect  So  1 ut  i  on  of  the  He  1  mho  1 tz  Equat  ion 

The  Helmholtz  equation  may  be  written  in  the  form 

delsqw-hw=f  <i> 

where  delsq  is  the  two-dimensional  Laplace  operator,  h  is  a  func- 
tion  of  y,   f  is  a  function  of  x  and  y,   and  w  is  the   dependent 
variable.    For  the  Finite  Element  (FE)  discretization,   consider 
the  grid  of  Fig.  1.   There  are  e  east-west  grid  lines  and  n  north- 
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Fig.  1.   Node  numbering  and  spacing. 

south  grid  lines,  defining  ne  intersections  (nodes).  Node  num- 
bering is  from  west  to  east,  beginning  at  the  southwest  corner. 
The  positive  x  direction  is  eastward  and  the  positive  y  direction 
is  northward.  In  the  Galerkin  FE  discretization  process  the  pro- 
duct h  w  is  treated  as  a  single  entity  and  the  tensor  product 
concept  is  employed  in  writing  the  matrix  form  of  the  equation. 
The  resu It  is 

MX  W  SY   +   SX  W  MY   -   MX  W  H  MY   =   V  <2> 

where  MX  and  SX  are  n  x  n  symmetric  "mass"  and  "stiffness"  mat- 
rices for  the  x  direction,  MY  and  SY  are  the  corresponding  e  x  e 
matrices  for  the  y  direction,  W  is  an  n  x  e  matrix  of  nodal 
values  of  the  dependent  variable,  H  is  diagonal  e  x  e,  and  V  is 
n  x  e.  MX  and  SX  depend  only  on  the  node  spacing  a(  and  MY  and 
SY  depend  only  on  the  node  spacing  bt  .  Defining  equations  for 
these  matrices  are  given  in  Appendix  A.  Note  that  matrices  SX 
and  SY  are  the  negatives  of  standard  FE  stiffness  matrices.    The 


columns  of  W  contain  nodal  values  of  w  from  the  west-east  rows, 
beginning  with  the  most  southerly.  If  the  nodal  values  of  f  are 
similarly  arranged  in  an  n  x  e  matrix  F,  then  V  is  the  product 
MX  F  MY.  The  nonzero  entries  in  the  diagonal  matrix  H  are  the 
values  of  h  for  the  successive  rows  of  nodes,  beginning  with  the 
most  souther  1 y . 

It   is   advantageous   for  the   succeeding   manipulations   to 
transpose  the  terms  of  <2>  to  obtain 

SY  WT  MX   +   MY  -WT  SX   -   MY  H  WT  MX   =   VT  <3> 

where  WT  and  VT  are  the  respective  transposes  of  W  and  V.  Note 
that  all  of  the  other  matrices  are  symmetric. 

Consider  the  generalized  eigenvalue  problem 

SX  Pi  =   l»  MX  pi  <4> 

where  pt  is  the  ith  eigenvector  (n  x  1)  and  lt  is  the  correspond- 
ing eigenvalue.  Assembling  the  eigenvectors  in  an  n  x  n  modal 
matrix  P  and  the  eigenvalues  in  an  n  x  n  (diagonal)  spectral 
matrix  L,  the  result  may  be  exhibited  as 

SX  P   =   MX  P  L  <5> 

It  is  advantageous  to  require  that  the  eigenvectors  be  normalized 
so  that  PT  MX  P  =  I,  where  PT  is  the  transpose  of  P  and  I  is  the 
nth  order  identity  matrix.  If  <3>  is  pos tmu 1 t i p 1 i ed  by  P  and  <5> 
is  used  to  replace  SX  P  in  the  second  terra,   the  result  is 

SY  WT  MX  P   +   MY  WT  MX  P  L   -   MY  H  WT  MX  P   =   VT  P     <6> 
Let  Q  =  WT  MX  P  and  substitute  in  <6>  to  get 

SY'  Q   +   MY  Q  L   =   U  <7> 

where  SY'  =  SY  -  MY  H  and  U  =  VT  P .  Observe  that  SY'  is  not 
symmetric,  but  is  tridiagonal.  Now  the  ith  column  qt  of  Q  may  be 
found  by  solving 


(SY'   +   1,  MY)  qt    =   u,         <8> 
where  ut  is  the  ith  column  of  U.    The  solutions  for  the  n  qt  may 
be  assembled  into  Q.   It  is  easy  to  show  that 

W   =   P  QT  <9> 

where  QT  is  the  transpose  of  Q.   This  completes  the  solution. 

Genera  1 i  zed  Ei  genva 1 ue  Prob 1  em  for  Symmetr  ic  Tr  idiagona 1  Matr  ices 
Solution  of  the  generalized  eigenvalue  problem  <'4>  for  the 
symmetric  tridiagonal  •  matrices  SX  and  MX  is  based  on  the  Sturm 
sequence  property  and  the  method  of  bisection.  The  theory  is 
given  by  Wilkinson  in  Ref.  6.  Using  this  theory  and  generalizing 
an  ALGOL  program  constructed  for  the  standard  eigenvalue  problem 
(Ref.  7),  SUBROUTINE  TR I E I G  has  been  written  to  find  the  eigen- 
va 1 ues . 

For  "any  eigenvalue,   the  corresponding  eigenvector   can   be 
found  by  rewriting  <4>  in  the  form 

(SX   -   14  MX)  qt   =   0  <10> 

If  the  first  component  of  qt  is  chosen  as  unity,  the  first  scalar 
equation  of  <10>  determines  the  second  component  and  each  suc- 
ceeding scalar  equation  determines  an  additional  component.  The 
vector  thus  found  may  then  be  normalized.  SUBROUTINE  EIGVEC  per- 
forms these  operations. 

Once  eigenvectors  have  been  determined,  SUBROUTINE  RAYLEE 
may  be  used  to  refine  the  eigenvalues  by  evaluating  the  Rayleigh 
quotient  for  each  eigenvector.  Following  this  an  additional  call 
to  SUBROUTINE  EIGVEC  enhances  the  accuracy  of  the  eigenvectors. 
Experience  to  date  indicates  that  further  iterations  are  super- 
f 1 uous . 


Genera  1 ized  Ei  genva 1 ue  Prob I  em  with  Per  iodic  Boundary  Cond  i  t ions 

When  there  is  a  periodic  boundary  condition  in  the  x  direc- 
tion, the  tridiagonal  form  of  matrices  SX  and  MX  is  modified  by 
the  presence  of  nonzero  elements  in  the  upper  right-hand  and 
lower  left-hand  corners.  The  matrices  remain  symmetric,  but  the 
Sturm  sequence  and  bisection  procedure  described  above  requires 
modification.  Although  the  underlying  strategy  is  unchanged, 
different  software  is  required.  In  this  instance,  the  eigen- 
values are  found  by  calling  SUBROUTINE  PEREIG.  To  determine  the 
eigenvectors  SUBROUTINE  EIGVCP  is  called.  A  new  problem  arises 
here,  because,  if  the  grid  spacing  is  uniform  in  the  x  direction, 
there  will  be  double  eigenvalues.  Specifically,  if  n,  the  number 
of  subdivisions  in  the  x  direction,  is  odd,  there  will  be  a 
single  zero  eigenvalue  and  all  of  the  others  will  be  double.  If 
n  is  even,  the  eigenvalue  of  largest  abso 1 ute* va 1 ue  will  also  be 
single.  Corresponding  to  double  eigenvalues,  the  eigenvectors 
are  not  unique  and  special  procedures  must  be  used  to  construct 
vector  pairs  having  the  required  orthogonality.  SUBROUTINE  EIG- 
VCP tests  for  repeated  eigenvalues  and  constructs  the  appropriate 
orthogonal  pairs  of  vectors  as  required.  For  this  periodic 
boundary  condition,  eigenvalue  refinement  via  the  Rayleigh  quo- 
tient is  effected  by  calling  SUBROUTINE  RAYLYP.  An  additional 
call  to  SUBROUTINE  EIGVCP  provides  better  eigenvectors.  Further 
iteration  is  believed  superfluous. 


Remarks 

Two  additional  subroutines  are  needed  to  utilize  the  pro- 
posed alterations  to  the  Stani f or th-Mi tche 1 1  program.  First  of 
these  is  SUBROUTINE  SETUP,  a  replacement  for  SUBROUTINE  EBVSET. 
The  second  is  SUBROUTINE  SLVBVP,  a  replacement  for  SUBROUTINE 
EBVP2D.   All  of  the  new  subroutines  are  listed  in  Appendix  A. 

The  routines  described  herein  have  been  tested  and  found  to 
perform  satisfactorily  on  a  12  x  12  grid.  The  parameter  RF  in 
TRIEIG  and  PEREIG  (currently  set  at  1.0  E-07)  may  need  to  be 
adjusted  for  larger  grids.  In  SUBROUTINE  EIGVCP  there  is  an 
additional  parameter,  SEP,  which  also  may  need  adjustment  for 
larger  grids.   It  is  currently  set  at  4.0  E-04. 

There  is  a  fundamental  problem  arising  from  the  use  of  single 
precision  arithmetic  in  solutions  to  the  Helmholtz  equation. 
This  results  from  the  fact  that  the  mean  geopotential  height  is 
approximately  three  orders  of  magnitude  greater  than  the  ampli- 
tude of  a  representative  disturbance.  Comparisons  between  single 
and  double  precision  solutions  show  clearly  that  round-off  ad- 
versely affects  the  former. 

Operation  counts  have  been  conducted  to  evaluate  the  potent- 
ial saving  resulting  from  substituting  the  Helmholtz  direct  sol- 
ver for  the  Concus-Golub  iterative  scheme.  For  a  13  x  13  grid 
and  assuming  only  one  cycle  of  iteration,  there  is  a  saving  of  10 
percent  for  each  time  step.  For  a  90  x  90  grid,  the  saving  is 
reduced  to  5  percent.  If  additional  iterations  are  required  for 
the  Concus-Golub  solution,  these  savings  would,  of  course,  be 
greater.  The  savings  resulting  from  the  proposed  eigenvalue  sol- 
utions have  not  been  evaluated,  but  are  expected  to  be  substan- 
tial . 
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APPENDIX  A.   "MASS"  and  MST I FFNESS"  MATRICES 
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Notes 


1.  Matrices  SX  and  SY  are  negatives  of  those  called  stiff 
ness  matrices  in  standard  FE  usage. 

2.  Forms  given  for  MX  and  SX  are  for  a  periodic  east-west 
boundary  condition.  For  a  Neumann  boundary  condition, 
the  terms  containing  a,  would  be  omitted. 


8 


APPENDIX  B.   FORTRAN  LISTINGS 


SUBROUTINE  SETUP ( E I GVEC, E I GVAL, B I GE, B I GC , B I GA , S , HX, HY , HELM , D I R , 
1  FOURTH, NI ,NJ, CON) 


*# 

* 

* 

*# 

C 

C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c- 


***************** 
SUBROUTINE  SETUP 
BOUNDARY  CONDITI 


***************************************************** 


IS  DESIGNED  TO  REPLACE  EBVSET  WHEN  THE  PERIODIC 
ON  IS  IMPOSED  ON  EASTERN  AND  WESTERN  BOUNDARIES 


********************************************************************** 
AUTHOR:  R.  E.  NEWTON,  SUMMER  1986 


ARGUMENTS 

OUT  -  EI GVEC 
EI GVAL 
BIGE 

BIGC 

BIGA 

IN  -     S 
HX 
HY 
HELM 
DIR 


MATRIX  OF  EIGENVECTORS,  N 
VECTOR  OF  EIGENVALUES,  NI 
DIAGONAL  FACTOR  OF  COEFFI 
DIRECTION 

SUPERD I  AGONAL  OF  UPPER  UN 
COEFFICIENT  MATRIX  FOR  Y 
SUBDIAGONAL  OF  LOWER  UNIT 
COEFFICIENT  MATRIX  FOR  Y 
SQUARE  OF  MAP  FACTOR,  NI 
NODE  SPACING  IN  X  DIRECT  I 
NODE  SPACING  IN  Y  D I RECT I 


LOGICAL 
LOGICAL 


SWITCH 
SWITCH 


FOURTH  -   LOGICAL  SWITCH  - 


TRUE  FOR 
TRUE  FOR 
POISSON 
TRUE  FOR 
POISSON 


I  X  NI ,   IN  X  DIRECTION 
,  IN  NOND€SCENDING  ORDER 
CIENT  MATRIX  FOR  Y 

IT  TRIANGULAR  FACTOR  OF 
DIRECTION 

TRIANGULAR  FACTOR  OF 
DIRECTION 
X  NJ 
ON 
ON 

HELMHOLTZ 

DIRICHLET 
PROBLEM 

FOURTH  ORDER 
PROBLEM 


PROBLEM 
B.C.  ON 


SOLUTION  OF 


NI 
NJ 


X-DIMENSION 
Y-DIMENSION 


NOTE 


N  IS  MAX(NI,NJ) 


PARAMETER(N=13) 

REAL  AS(N) ,BS(N) ,CS(N) , S(NI ,NJ),HX(NI ) ,HY(NJ) , 

1  AR(N) ,BR(N) ,CR(N),BIGA(NI , NJ ) ,BIGC(NI , NJ ) ,BIGE(NI , NJ ) 

2  E I GVEC (N I , NI ) ,EIGVAL(NI ) , AM ( N ) , BM ( N ) , CM ( N ) , CON 
LOGICAL  HELM, DIR, FOURTH 


WU(N) 


CF  =  2. 
IF(FOURTH)  CF  = 


5. 


C 
C 

c 


c 
c 
c 


SET  UP  "MASS"  MATRIX  FOR  X  DIRECTION  (SUBDIAGONAL  AM,  DIAGONAL  BM, 
SUPERDIAGONAL  CM).  PERIODIC  BOUNDARY  CONDITION  ASSUMED.  IF  NOT 
PERIODIC,  CALL  SETABC  INSTEAD. 

CALL  SETABX(AM,BM,CM,HX,CF,NI ) 

SET  UP  "STIFFNESS"  MATRIX  FOR  X  DIRECTION  (SUBDIAGONAL  AS,  DIAGONAL 
BS,  SUPERDIAGONAL  CS).   PERIODIC  BOUNDARY  CONDITION  ASSUMED.    IF 
NOT  PERIODIC,  CALL  SETD2  INSTEAD. 

CALL  SETD2N(AS, BS,CS, 1. ,HX, NI ) 

SOLVE  THE  GENERALIZED  EIGENVALUE  PROBLEM:  K  X  =  Z  M  X,  WHERE  K  AND 


c 
c 
c 


c 
c 

c 


M  ARE  STIFFNESS  AND  MASS  MATRICES, 
VECTOR  CORRESPONDING  TO  EIGENVALUE 
ASSUMED.    IF  NOT  PERIODIC,  CALL  TR 


RESPECTIVELY 
Z.  PERIODIC 
EIG  INSTEAD. 


AND  X  IS  THE  EI  GEN- 
BOUNDARY  CONDITIONS 


CALL  PEREIGCEIGVEC, E I GVAL, AS , BS , BM, AM, WU , N I ) 


SET  UP  "MASS"  MATRIX  (SUBDIAGONAL  AM,  DIAGONAL  BM,  SUPERD I  AGONAL  CM! 
AND  "STIFFNESS"  MATRIX  (SUBDIAGONAL  AS,  DIAGONAL  BS,  SUPERD I  AGONAL 
CS)  FOR  THE  Y  DIRECTION. 

CALL  SETABC(AM,BM,CM,HY,CF,NJ) 
CALL  SETD2(AS,BS,CS, 1. ,HY,NJ) 
IF( .NOT.HELM)GO  TO  11 


C   FOR  THE  HELMHOLTZ  PROBLEM  THE  STIFFNESS  MATRIX  IS  MODIFIED  BY  SUB 
C   TRACT  I NG  THE  CORRESPONDING  ELEMENT  OF  THE  MASS  MATRIX  DIVIDED  BY 
C   (GZMN  TIMES  DELT**2  TIMES  THE  RELEVANT  SQUARE  OF  THE  MAP  FACTOR). 
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AS(1  ) 
BS(1) 
CS(1) 
NJM  = 
DO  10 

AS  (J) 
BS(J) 
CS(J  ) 
CONTINUE 
AS(NJ)  = 
BS(NJ)  = 
CS(NJ)  = 


=  0. 
=  BS(1 ) 
=  CS(1) 
NJ  -  1 
J=2, NJM 


-  BM(1 )*C0N/S(1, 1) 

-  CM( 1 )*CON/S( 1,2) 


=  AS(J)  -  AM( J )*CON/S(i, J-l ) 
=  BS(J)  -  BM( J)*C0N/S(1, J) 
=  CS(J)  -  CM(J )*C0N/S(1, J+l) 


AS(NJ ) 
BS(NJ) 
0. 


-  AM(NJ)*CON/S(i, NJM) 

-  BM(NJ)*CON/S( 1,NJ) 


C   LOOP  OVER  EIGENVALUES  TO  CONSTRUCT  B I GA ,  B I GC ,  AND  B I GE 


11 


12 


DO  20  I  =  1,  NI 

EIG  =  EI GVAL ( I ) 


DO 


12  J  =  1, 

AR(J)  = 

BR(J)  = 

CR(J)  = 
CONTINUE 
IF( .NOT. HELM. AND 


NJ 

AS(J  ) 
BS(  J) 
CS(J  ) 


EIG*AM< J) 
EIG*BM<J) 
EIG*CM(J) 


NOT.DIR. AND. I .EQ.NI ) CR ( 1 )  =  0 


C   FACTOR  COEFFICIENT  MATRIX 

CALL  SETTRI ( BR , CR, AR, AR, BR, CR , NJ ) 

C   STORE  RESULTS  IN  BIGE,  B I GC  AND  BIGA.   NOTE  SIGN  CHANGES  FOR  B I GC 
C   AND  BIGA. 


16 
20 


DO  16  J  =  1,  NJ 
BIGE( I , J) 
BIGC( I , J ) 
BIGA(  I  , J  ) 
CONTINUE 
CONTINUE 
RETURN 
END 


BR(  J) 
-CR(J ) 
-AR( J) 
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cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  PERE I G ( V, X, B, C, D, E, WU, N) 


a*********************************************************************** 

*  SUBROUTINE  PEREIG  USES  THE  METHOD  OF  BISECTION  TO  FIND  EIGENVALUES   * 

*  FOR  THE  GENERALIZED  EIGENVALUE  PROBLEM  INVOLVING  SYMMETRIC  (ALMOST)  * 

*  TRIDIAGONAL  MATRICES  WITH  NONZERO  ELEMENTS  IN  'CORNERS'  (PERIODIC    * 

*  BOUNDARY  CONDITION).  IT  IS  ADAPTED  FROM  AN  ALGOL  PROGRAM  NAMED       * 

*  'BISECT'  WRITTEN  BY  BARTH,  MARTIN,  AND  WILKINSON  (NUM.  MATH.  9,  386-* 

*  393  (1967)).   'BISECT'  FINDS  EIGENVALUES  (STANDARD  EIGENVALUE  PROB-  * 

*  LEM)  FOR  A  SYMMETRIC  TRIDIAGONAL  MATRIX.   THE  MATRICES  FOR  PEREIG    * 

*  ARE  INPUT  AS  VECTORS  -  C(N)  AS  THE  DIAGONAL  AND  B(N)  AS  THE  SUB-     * 

*  DIAGONAL  OF  THE  "STIFFNESS"  MATRIX  AND  D(N)  AS  THE  DIAGONAL  AND  E(N)* 

*  AS  THE  SUBDI  AGONAL  OF  THE  "MASS"  MATRIX.   B(l)  AND  Ed)  ARE  THE      * 

*  "CORNER"  ELEMENTS  OF  THE  RESPECTIVE  MATRICES.   SUBROUTINE  EIGVCP     * 

*  FINDS  THE  EIGENVECTORS  AND  NORMALIZES  THEM  WITH  RESPECT  TO  THE  MASS  * 

*  MATRIX.   SUBROUTINE  RAYLYP  USES  THR  RAYLEIGH  QUOTIENT  TO  OBTAIN  IM-  * 

*  PROVED  EIGENVALUES  USING  THESE  EIGENVECTORS.   A  SECOND  CALL  TO       * 

*  TO  EIGVCP  EFFECTS  A  CORRESPONDING  IMPROVEMENT  IN  THE  EIGENVECTORS.  * 
a*********************************************************************** 
C 

E.  NEWTON,  SUMMER  1986 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


AUTHOR:  R. 


ARGUMENTS 
OUT  - 


IN 


V 

X 

C 

B 

D 

E 

WU 

N 


MATRIX  OF  EIGENVECTORS,  N  X  N. 

WITH  RESPECT  TO  "MASS"  MATRIX 

VECTOR  OF  EIGENVALUES  IN  NONDESCEND ING 

DIAGONAL  OF  "STIFFNESS"  MATRIX 

OF  "STIFFNESS"  MATRIX 

"MASS"  MATRIX 

OF  "MASS"  MATRIX 


NORMALIZED  WITH 


ORDER 


SUBDI AGONAL 
DIAGONAL  OF 
SUBDIAGONAL 
WORK  VECTOR 
VECTOR  SIZE 


(=  NI  ) 


NOTE 


MATRICES  ARE  FOR 
CONDITION.   NODE 


X-DIRECTION  WITH 
SPACING  IS  HX. 


PERIODIC  BOUNDARY 


INTEGER  N,Z, I , A,K,N1 

REAL  C(N),B(N),X(N),WU(N) , EPS1 , EPS2, RF, Fl , XM I N, XMAX, 

1  Xl,XU,XO,D(N),E(N) , DC,DB,DD,DE,TMAX,TMIN, V(N,N) , 

2  Q,R,S,Q1,R1,R2, DR, QTEMP , ROLD, SW 
DATA  EPSl,RF/0. , l.E-7/ 

N1=0 
Z  =  0 


C   CALCULATION  OF  XMAX,  XM I N 

DC=C(N) 
DB=ABS(B(N) ) 
DD=D(N) 
DE=E(N) 

XMAX= (DC+DB) / (DD+DE) 

XMIN= (DC-DB)/ (DD-DE) 

C   EIGENVALUES  ASSUMED  NEGATIVE. 


IF  EIGENVALUES  ARE  ALL  POSITIVE, 
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C   REPLACE  THE  TWO  PRECEDING  LINES  WITH  THE  FOLLOWING  TWO  LINES. 
»      XMAX=(DC+DB) /(DD-DE) 

*  XMIN=(DC-DB)/(DD+DE) 
DO  2  I-N-l, 1, -1 

DC  =  C(  I  ) 

DB  =  ABS(B(  I  )  )+ABS(B(  I  +1  )  ) 

DD=D( I ) 

DE=E( I )+E( 1+1 ) 

TMAX=(DC+DB)/(DD+DE) 

TMIN= (DC-DB)/ (DD-DE) 
C   EIGENVALUES  ASSUMED  NEGATIVE.    IF  EIGENVALUES  ARE  ALL  POSITIVE, 
C   REPLACE  THE  PRECEDING  TWO  LINES  WITH  THE  FOLLOWING  TWO  LINES. 

*  TMAX=(DC+DB)/(DD-DE) 

*  TMIN=(DC-DB)/(DD+DE) 

I F ( TMAX . GT . XMAX ) XMAX  =  TMAX 
I F ( TM I N . LT . XM I N ) XM I N  =  TM I N 
2      CONTINUE 

C   SET  EPS2 

IFCXMIN+XMAX.GT.O. ) THEN 

EPS2=RF*XMAX 
ELSE 

EPS2=RF*(-XMIN) 
END  IF 

IF(EPSl.LE.O. )EPS1=EPS2 
EPS2=0.5*EPSl+7. *EPS2 

C   INNER  BLOCK 

XO=XMAX 
DO  4  1=1, N 

X( I )=XMAX 

WU( I )=XMIN 
4      CONTINUE 

C   LOOP  FOR  K-TH  EIGENVALUE 

DO  100  K=N, 1, -1 
XU=XMIN 
DO  6  I=K, 1, -1 

IF(XU.LT. WU( I ) )THEN 
XU  =  WU(  I  ) 
GO  TO  10 
END  IF 
6  CONTINUE 

10         IFCXO.GT.X(K) )XO=X(K) 
20         X1=(XU+X0)/2.D0 
Z  =  Z  +  1 

C   STURM'S  SEQUENCE 

A  =  0 

Q=C(1)-X1*D(1) 
Q_1=C(2)-X1*D<2) 
R=B( 1 ) -X1*E< 1 ) 
S=C(N)-X1*D(N) 
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R1=B(2)-X1*E(2) 

R2  =  0. 

DR  =  0. 

DO  30  1=2, N-2 

IF(Q.EQ.O. )THEN 
SW=1. 

IF(Q1+2.*R1.EQ.0.)SW=-1. 
Q=Q1+SW*2. *R1 
R1=R1+SW*Q1 

R2=SW*(B(  I+1)-X1*E(I+D) 
END  IF 

IFCQ.LT.O. )A=A+1 
S=S-R*R/Q 
ROLD=R 
R=DR-R1*R/Q 
QTEMP  =  Q_1-R1*R1/Q 
R1=B(I+1)-X1*E( I+1)-R1*R2/Q 
Qi=C(I+l)'-Xl*D(  I+1)-R2*R2/Q 
DR=-R0LD*R2/Q 
R2=0.D0 
Q=QTEMP 
30       CONTINUE 

IF(Q.EQ.O.DO)  THEN 
SW=1.D0 

JF<Q1+2.D0*R1.EQ.0.D0)SW=-1.D0 
Q=Q1+2.D0*SW*R1 
R1=R1+SW*Q1 
R=R+SW*(B(N)-X1*E(N) ) 
IF(Q.LT.O.DO)A=A+i 
ELSE 

IFCQ.LT.O. )A=A+1 
END  IF 
S=S-R*R/Q 

R=B(N)-X1*E(N)-R*R1/Q 
Q1=Q1-R1*R1/Q 
IFCQ1.EQ.0. AND.R.NE.O. )  THEN 

A  =  A+1 
ELSE 

IF(Ql.LT.O. )A=A+1 
IF(S-R*R/Ql.LT.O. )A=A+1 
END  IF 
IF(A.LT.K)THEN 

IF(A.LT.  1)THEN 
WU(1)=X1 
XU  =  X1 
ELSE 

WU(A+1)=X1 
XU  =  X1 

IF(X(A) .GT.X1)X(A)=X1 
END  IF 
ELSE 

X0  =  X1 
END  IF 

IF(X0-XU.GT.2. *RF*(ABS(XU)+ABS(XO) )+EPSi)GO  TO  20 
X(K)=(X0+XU)/2. 
100    CONTINUE 

CALL  EIGVCP(V,B,C,D,E,X,N) 
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CALL  RAYLYPCX, V,B,C, WU , N ) 
CALL  EIGVCP(V,B,C,D, E,X,N) 
RETURN 
END 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx> 

SUBROUTINE  E I GVCP ( V, B, C, D, E, X, N ) 


ftK#ft#ft«»»*######ft*****************************************************» 

*  SUBROUTINE  EIGVCP  FINDS  EIGENVECTORS  BY  DIRECT  SOLUTION  OF  THE  GOV- 

*  ERNING  LINEAR  ALGEBRAIC  EQUATIONS.   FOR  ANY  PAIR  OF  EQUAL  EIGEN- 

*  VALUES,  THE  CORRESPONDING  VECTORS  ARE  MADE  ORTHOGONAL  WITH  RESPECT 

*  TO  THE  "MASS"  MATRIX.   ALL  VECTORS  ARE  NORMALIZED  WITH  RESPECT  TO 

*  THE    "MASS"    MATRIX. 

***********************  N  ********  M  ***************  ft  ft  ft  ft  ft  ft  ft  ft  ***************  4 
C 

AUTHOR   -   R.  E.  NEWTON,  SUMMER  1986 


ARGUMENTS 
OUT  - 
IN   - 


V 
B 
C 
D 
E 
X 
N 


MODAL  MATRIX, 
SUBD I  AGONAL 
DIAGONAL  OF 
DIAGONAL  OF 
SUBDIAGONAL 
VECTOR  OF  E 
VECTOR  SIZE 


N  X  N 


(COLUMNS  ARE  EIGENVECTORS 


OF  STIFFNESS  MATRIX 
STIFFNESS  MATRIX 
MASS  MATRIX 
OF  MASS  MATRIX 
GENVALUES 
(=NI  ) 


INTEGER  J,K,N,L,N1 
PARAMETER(N1=12) 

REAL  V(N,N),P(N1),Q(N1),T(N1),X(N),B(N),C(N),D(N),E(N), 
1      H, VN,TN,DIAG,X1,X2, VTMV , TTMV , TTMT, SUM 


L  =  0 
DO  30 
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K=1,N 
IF(L.NE.O)THEN 

L  =  0 

GO  TO  30 
END  IF 
X1=X(K) 
V( 1,K)=1. DO 
V(2,K)=0.D0 
T(1)=0.D0 
T(2)=1.D0 
P(  D=B(l)-XlftE(l) 
P(2)=B(2)-XlftE(2) 
VN=-(C( 1)-X1*D(1))/P(1) 
TN=-P(2)/P(1) 
DO  10  J=2,N-1 

DIAG=C(J)-X1«D(J) 

P(J+1)=B(J+1)-X1*E(J+1) 

V(J+l,K)=-(P(J)*V(J-l,K)+DIAG*V(J,K))/P(J+i) 

T(J+1)=-(P(J)*T(J-1)+DIAG»T(J))/P(J+1) 
CONTINUE 
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C   CHECK  FOR  A  DOUBLE  ROOT 

IF(K.EQ.N)GO  TO  101 

X2=X(K+1) 

CRIT  =  ABS(  (X2-XD/CX2  +  X1)  ) 

IF(CRIT.LT.4.E-4)G0    TO    22 

C   CONSTRUCT  EIGENVECTOR  FOR  SINGLE  ROOT 

101       H=-(V(N,K)-VN)/(T(N)-TN) 
DO  11  J=2,N 

V(J,K)=V(J,K)+H*T(J) 

11  CONTINUE 

C   NORMALIZE  WITH  RESPECT  TO  MASS  MATRIX 

P(1)=D(1)*V(1,K)+E(2)*V(2,K)+E(1)*V(N,K) 
DO  12  J=2,N-1 

P< J)=E(J)*V(J-1,K)+D(J)*V(J,K)+E(J+1)*V(J+1,K) 

12  CONTINUE 
P(N)=E(N)*V(N-l,K)+D(N)*V(N,K)+E(l)*V(i,K) 
VTMV=O.DO 

DO  14  J=1,N 

VTMV=VTMV+P(J)*V(J,K) 
14        CONTINUE 

VTMV= 1 . DO/SQRT ( VTMV ) 
DO  16  J=1,N 

V(J,K)=V(J,K)*VTMV 
16        CONTINUE 
GO  TO  30 

C   CONSTRUCT  EIGENVECTORS  FOR  DOUBLE  ROOT  AND  NORMALIZE 

22  L=l 
P(1)=D(1)*V(1,K)+E(2)*V(2,K)+E(1)*V(N,K) 
Q(l)=D(l)*T(i)+E(2)*T(2)+E(l)*T(N) 

DO  23  J=2,N-1 

P<J)=E(J)*V(J-l,K)+D<J)*V(J,K)+E(J+i)*V(J+l,K) 
Q(J)=E(J)*T(J-1)+D<J)*T(J)+E<J+1)*T(J+1) 

23  CONTINUE 
P(N)=E(N)*V(N-1,K)+D(N)*V(N,K)+E(1)*V(1,K) 
Q(N)=E(N)*T(N-1)+D(N)*T(N)+E(1)*T(1) 
VTMV=0. DO 

TTMV=O.DO 
TTMT=O.DO 
DO  24  J=1,N 

VTMV=VTMV+P( J)*V( J,K) 

TTMV=TTMV+P( J)*T(J) 

TTMT  =  TTMT  +  Q( J  >*T< J ) 

24  CONTINUE 
H=-TTMV/VTMV 

SUM=TTMT+2. DO*H*TTMV+H*H*VTMV 

VTMV =1. DO/SQRT (VTMV) 

SUM=1. DO/SQRT(SUM) 

H=H*SUM 

DO  26  J=1,N 

V(J,K+1)=V(J,K)*H+T(J) *SUM 
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V(J,K)=V(J,K)*VTMV 
26         CONTINUE 
30     CONTINUE 

RETURN 

END 


cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  RAYLYP ( X, V, B, C, P, N ) 


*  # 
** 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c- 


SUBROUTINE  RAYLYP  USES  THE  RAYLEIGH  QUOTIENT  TO  FIND  IMPROVED 
EIGENVALUES  FROM  THE  ALREADY  NORMALIZED  EIGENVECTORS 

a******************************************************************** 


AUTHOR 

ARGUMENTS 
OUT  - 
IN   - 


R.  E.  NEWTON,  SUMMER  1986 


X  -  VECTOR  OF  N  EIGENVALUES  IN  NONDESCEND I NG  ORDER 

V  -  MODAL  MATRIX,  N  X  N.  (NORMALIZED  WITH  RESPECT  TO 

MASS  MATRIX. ) 

B  -  SUBD I  AGONAL  OF  STIFFNESS  MATRIX 

C  -  DIAGONAL  OF  STIFFNESS  MATRIX 

P  -  WORK  VECTOR 

N  -  VECTOR  SIZE  (=  NI ) 


INTEGER  J,K,N 

REAL  V(N,N) ,P(N) ,X(N) ,B(N) ,C(N) ,X1 

DO  20  K=1,N 

P(1)=C(1)*V(1,K)+B(2)*V(2,K)+B(1)*V(N,K) 
DO  12  J=2,N-1 

P(J)=B(J)*V(J-1,K)+C(J)*V(J,K)+B(J+1)*V(J+1,K) 
12         CONTINUE 

P(N)=B(N)*V(N-1,K)+C(N)*V(N,K)+B(1)*V(1,K) 

X1=0.D0 

DO  16  J=1,N 

X1=X1+P(J)*V(J,K) 
16         CONTINUE 

X(K)=X1 
20     CONTINUE 
RETURN 
END 

■ 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  SLVBVPCPHI , RHS , E I GVEC , B I GE, B I GC, B I GA, WK, DIR,HELM, 
1  NI ,NJ ) 

a********************************************************************* 

*   SUBROUTINE  IS  A  SUBSTITUTE  FOR  EBVP2D 

********************************************************************** 

C 

C   AUTHOR      -     R.  E.  NEWTON,  SUMMER  1986 

C 

C   ARGUMENTS 
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c 

OUT  - 

PHI 

c 

IN   ■ 

RHS   - 

c 

EIGVEC  - 

c 

BIGE   - 

c 

BIGC   - 

c 

c 

BIGA   - 

c 

c 

WK 

c 

DIR   - 

c 

HELM   - 

c 

NI 

c 
c--- 

NJ 

SOLUTION 

RIGHT-HAND  SIDE 

MATRIX  OF  EIGENVECTORS 

INVERSE  OF  DIAGONAL  FACTOR  OF  COEFFICIENT  MATRIX 

SUPERD I  AGONAL  OF  UNIT  UPPER  TRIANGULAR  FACTOR  OF 

COEFFICIENT  MATRIX 

SUBD I  AGONAL  OF  UNIT  LOWER  TRIANGULAR  FACTOR  OF 

COEFFICIENT  MATRIX 

WORK  AREA 

LOGICAL  SWITCH  -  TRUE  FOR  DIRICHLET  B.C. 

LOGICAL  SWITCH  -  TRUE  FOR  HELMHOLTZ  EQUATION 

X-DIMENSION 

Y-DIMENSION 


REAL  PHI (NI ,NJ) ,RHS(NI , NJ ) , EIGVEC (N I , N I ) , B I GE(N I , NJ ) , B I GC (N I , NJ ) , 
1   BIGACNI ,NJ) , WK(NJ,NI ) , DUM 
LOGICAL  DIR,  HELM 

C   TRANSFORM  RIGHT-HAND  SIDE 

CALL  MATPR(WK,RHS,EIGVEC,NJ,NI , N I , N J , N I , N I , DUM, DUM, DUM, D I R, 
1  .FALSE. ) 

C   PERFORM  FORWARD  REDUCTIONS 

DO  5  I  =  1,  NI 

WK(1, I )  =  WK(i, I )*BIGE( 1,1) 
5      CONTINUE 

.DO  15  I  =  1,  NI  . 
DO  10  J  =  2,  NJ 

WK(J,I)  =  WK(J, I )*BIGE( I , J)  +  WK( J-l, I )*BIGA( I , J) 
10        CONTINUE 
15     CONTINUE 

C   PERFORM  BACK  SUBSTITUTIONS 

DO  25  I  =  1,  NI 

DO  20  J  =  NJ-1,  1,  -1 

WK(J,I)  =  WK(J,I)  +  WK(J+1, I )*BIGC( I , J) 
20        CONTINUE 
25     CONTINUE 

C   BACK  TRANSFORM  RESULTS 

DO  40  J  =  1,  N J 

DO  35  I  =  1,  NI 
DUM  =  0. 
DO  30  K  =  1,  NI 

DUM  =  DUM  +  EIGVEC( I ,K)*WK(J,K) 
30  CONTINUE 

PHI ( I , J)  =  DUM 
35        CONTINUE 
40     CONTINUE 
RETURN 
END 
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cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  TRIEIG(V,X,B,C,D,E,WU,N) 


*  SUBROUTINE  TRIEIG  USES  THE  METHOD  OF  BISECTION  TO  FIND  EIGENVALUES  » 

*  FOR  THE  GENERALIZED  EIGENVALUE  PROBLEM  INVOLVING  SYMMETRIC  TRIDIA- 

*  GONAL  MATRICES.   THE  ROUTINE  IS  ADAPTED  FROM  AN  ALGOL  PROGRAM  NAMED* 

*  'BISECT'  WRITTEN  BY  BARTH,  MARTIN,  AND  WILKINSON  (NUM.  MATH.  9,  386- 

*  393  (1967)).    'BISECT'  FINDS  EIGENVALUES  (STANDARD  EIGENVALUE  PROB- 

*  LEM)  FOR  A  SYMMETRIC  TR I D I  AGONAL  MATRIX.   THE  MATRICES  FOR  TRIEIG 

*  ARE  INPUT  AS  VECTORS  -  C(N)  AS  THE  DIAGONAL  AND  B(N)  AS  THE  SUB-    * 

*  DIAGONAL  OF  THE  "STIFFNESS"  MATRIX  AND  D(N)  AS  THE  DIAGONAL  AND  E(N) 

*  AS  THE  SUBDIAGONAL  ON  THE  "MASS"  MATRIX.   SUBROUTINE  EIGVEC  FINDS 

*  THE  EIGENVECTORS  AND  NORMALIZES  THEM  WITH  RESPECT  TO  THE  MASS  MATRIX 

C 

AUTHOR:  R.  E.  NEWTON,  SUMMER  1985 


C 

C 

c 
c 

c 
c 
c 
c 
c 
c 
c 
c 
c 


ARGUMENTS 
OUT  - 


IN  - 


X 

C 

B 

D 

E 

WU 

N 


MATRIX  OF  EIGENVECTORS,  N  X  N, 
RESPECT  TO  THE  "MASS"  MATRIX 
VECTOR  OF  EIGENVALUES  IN  NONDESCEND I NG 
DIAGONAL  OF  "STIFFNESS"  MATRIX 

OF  "STIFFNESS"  MATRIX 

"MASS"  MATRIX 

OF  "MASS"  MATRIX 


NORMALIZED  WITH 


ORDER 


SUBDIAGONAL 
DIAGONAL  OF 
SUBDIAGONAL 
WORK  VECTOR 
VECTOR  SIZE 


(=  NI  ) 


INTEGER  N,M1,N,Z, I , A,K 

REAL  C(N),B(N),X(N),WU(N) , EPS1 , EPS2, RF, Fl , XM I N , XMAX , 
1Q,X1,XU,X0,DELT,D(N) ,E(N) , DC, DB , DD, DE, TMAX , TM I N , V ( N , N ) ,G(N,N) 
DATA  EPSl,RF/0. , l.E-7/ 

C   CALCULATION  OF  XMAX,  XM I N 

B( 1)=0. 

E(1)=0. 

DC=C(N) 

DB=ABS(B(N) ) 

DD=D(N) 

DE=E(N) 

XMAX=(DC+DB)/( DD+DE) 

XMIN= (DC-DB)/ (DD-DE) 
C   NEGATIVE  EIGENVALUES  ASSUMED.    IF  EIGENVALUES  ARE  ALL  POSITIVE, 
C   REPLACE  THE  TWO  PRECEDING  LINES  WITH  THE  FOLLOWING  TWO  LINES. 

*  XMAX=(DC+DB) / (DD-DE) 

*  XMIN=(DC-DB) / (DD+DE) 
DO  2  I=N-1, 1, -1 

DC=C( I ) 

DB  =  ABS(B(  I))+ABS(B(I+D) 

DD=D( I ) 

DE=E( I ) +E( I +1 ) 

TMAX=(DC+DB) / (DD+DE) 
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TMIN=(DC-DB)/(DD-DE) 
C   NEGATIVE  "EIGENVALUES  ASSUMED.    IF  EIGENVALUES  ARE  ALL  POSITIVE, 
C   REPLACE  THE  TWO  PRECEDING  LINES  WITH  THE  FOLLOWING  TWO  LINES. 

*  TMAX=(DC+DB)/ (DD-DE) 

*  TMIN=(DC-DB)/(DD+DE) 
IF(TMAX.GT.XMAX)XMAX=TMAX 
IF(TMIN.LT.XMIN)XMIN=TMIN 

2      CONTINUE 

C   SET  EPS2 

IFCXMIN+XMAX.GT.O. ) THEN 

EPS2=RF*XMAX 
ELSE 

EPS2=RF*(-XMIN) 
END  IF 

IF(EPSl.LE.O. )EPS1=EPS2 
EPS2=Q.5*EPSl+7. *EPS2 

C   INNER  BLOCK 

XO=XMAX 

DO  4  I=N, 1, -1 

X( I )=XMAX 

WU( I )=XMIN 
4      CONTINUE 

C   LOOP  FOR  K-TH  EIGENVALUE 

DO  100  K=N, 1,-1   • 
XU=XMIN 
DO  6  I=K, 1, -1 

IF(XU.LT. WU( I ) )THEN 
XU=WU( I ) 
GO  TO  10 
END  IF 
6      CONTINUE 
10        IF(XO.GT.X(K) )XO=X(K) 
20        Xl=(XU+X0)/2. 
Z  =  Z+1 

C   STURM'S  SEQUENCE 

A  =  0 

Q=l. 

DO  30  I=1,N 

IF(Q.EQ.O. )THEN 

DELT=ABS(B( I ) -XI *E ( I ) ) /RF 
ELSE 

DELT=(B( I )-Xl*E( I ) )**2/Q 
END  IF 

Q  =  C(  I  )-Xl*D< I )-DELT 
IFCQ.LT.O. )A=A+1 
30       CONTINUE 

IF(A.LT.K)THEN 

IF(A.LT. 1 )THEN 
WU(1)=X1 
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XU  =  X1 
ELSE 

WU(A+i )=X1 
XU  =  X1 

IF(X(A) .GT.X1 )X(A)=X1 
END  IF 
ELSE 

X0  =  X1 
END  IF 
IF(XO-XU.GT. 2. *RF*(ABS(XU) +ABSCXO) ) +EPS1 )G0  TO  20 
X(K)=(X0+XU)/2. 
100    CONTINUE 

CALL  EIGVEC(V,B,C,D,E,X,N) 
CALL  RAYLEE(V,B,C,X,N) 
CALL  EIGVEC(V,B,C,D,E,X,N) 
RETURN 
END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX) 

SUBROUTINE  E I GVEC ( V, B, C, D, E, X, N ) 

*  SUBROUTINE  E I GVEC  FINDS  EIGENVECTORS  BY  DIRECT  SOLUTION  OF  THE  GOV- 

*  ERNING  LINEAR  ALGEBRAIC  EQUATIONS  AND  NORMALIZES  THEM  WITH  RESPECT 

*  TO  THE  MASS  MATRIX, 
a********************************************************************** 

C 

C   AUTHOR      -   R.  E.  NEWTON,  SUMMER  1985 

C 

C   ARGUMENTS 

C      OUT      V         MODAL  MATRIX,  N  X  N.    (COLUMNS  ARE  EIGENVECTORS.: 

C       IN   -    B     -   SUBD I  AGONAL  OF  STIFFNESS  MATRIX 

C  C     -   DIAGONAL  OF  STIFFNESS  MATRIX 

C  D     -   DIAGONAL  OF  MASS  MATRIX 

C  E     -   SUBDIAGONAL  OF  MASS  MATRIX 

C  N         VECTOR  SIZE  (=  NI ) 

C 


INTEGER  J,K,N 

REAL  V(N,N),P(5),X(N),B(N),C(N),D(N),E(N) ,X1,DSQRT, SUM 

DO  20  K=1,N 
X1=X(K) 
V(1,K)=1. 
P(2)=B(2)-E(2)*X1 
V(2,K)=(D(1)*X1-C(1))/P<2) 
DO  10  J=2,N-1 

P(J+1)=B(J+1)-E(J+1)*X1 

V(J+1,K)=-(P(J)*V(J-1,K)+(C(J)-D(J)*X1)*V(J,K))/P(J+1) 
10         CONTINUE 

C   NORMALIZE  WITH  RESPECT  TO  MASS  MATRIX 

P(1)=D(1)*V(1,K)+E(2)*V(2,K) 
DO  12  J=2,N-1 


20 


12 


14 


16 
20 


P(J)=E(J)*V(J-1,K)+D(J)*V(J,K)+E(J+1)*V(J+1,K) 
CONTINUE 

P(N)=E(N)*V(N-1,K)+D(N)*V(N,K) 
SUM=0. 
DO  14  J=1,N 

SUM=SUM+P(J)*V(J,K) 
CONTINUE 
SUM=1. /SQRT(SUM) 
DO  16  J=1,N 

V( J,K)=V(J, K)*SUM 
CONTINUE 
CONTINUE 
RETURN 
END 


cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  RAYLEE-<  V,  B,  C,  X,  N> 


************************************************************************ 

*  SUBROUTINE  RAYLEE  USES  THE  RAYLIEGH  QUOTIENT  TO  FIND  IMPROVED  EIGEN-* 

*  VALUES  FROM  THE  ALREADY  NORMALIZED  EIGENVECTORS.  * 
************************************************************************ 

C 

AUTHOR      -   R.  E.  NEWTON,  SUMMER  1985 


C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 


ARGUMENTS 
OUT  - 
IN   - 


X 

V 

B 
C 
N 


VECTOR  OF  EIGENVALUES  IN  ASCENDING  ORDER 

MODAL  MATRIX,  N  X  N.    (NORMALIZED  WITH  RESPECT  TO 

MASS  MATRIX) 

SUBD I  AGONAL  OF  STIFFNESS  MATRIX 

DIAGONAL  OF  STIFFNESS  MATRIX 

VECTOR  SIZE  (=  NI ) 


INTEGER  J,K,N 

REAL  V(N,N),P(5),X(N),B(N),C(N),X1 

DO  20  K=1,N 

P(1)=C(1)*V(1,K)+B(2)*V(2,K) 
DO  12  J=2,N-1 

P(J)=B(J)*V(J-1,K)+C(J)*V(J,K)+B(J+1)*V(J+1,K) 
12        CONTINUE 

P(N)=B(N)*V(N-1,K)+C(N)*V(N,K) 

X1=0. 

DO  16  J=1,N 

X1=X1+P( J) *V( J,K) 
16        CONTINUE 

X(K)=X1 
20     CONTINUE 
RETURN 
END 
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